ELSEVIER 


Available  online  at  www.sciencedirect.com 


SCIENCE 


4 


DIRECT® 


Journal  of  Power  Sources  118  (2003)  79-85 


JOURNAL  OF 


www.elsevier.com /locate /jpowsour 


Computer  methods  for  performance  prediction  in  fuel  cells 

S.B.  Bealea’*,  Y.  Lina,  S.V.  Zhubrinb,  W.  Donga  l 

3 National  Research  Council,  Montreal  Road  Ottawa,  Ottawa,  Ont.,  Canada  K1A  0R6 
b Concentration  Heat  and  Momentum  Ltd.,  40  High  Street,  Wimbledon  Village,  London  SW19  5AU,  UK 


Abstract 

Several  transport  models  for  fuel  cells  have  been  developed.  The  models  are  compared  and  tested  for  single  fuel  cells  and  multi-cell  stacks 
of  planar  solid-oxide  fuel  cells,  the  three  main  approaches  considered  are  (a)  a  detailed  numerical  model  (DNM)  of  flow,  heat  and  mass 
transfer  and  electrochemistry,  (b)  a  flow-based  methodology  based  on  a  distributed  resistance  analogy  (DRA),  and  (c)  a  presumed-flow 
methodology  (PFM).  The  results  from  each  of  the  above  approaches  are  compared  in  detail,  and  merits  and  drawbacks  discussed.  It  is  shown 
that,  under  certain  circumstances,  the  simpler  approaches  have  the  potential  to  supplant  or  complement  the  direct  numerical  method  in  the 
analysis  of  fuel  cells. 
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1.  Introduction 

High  power  generation  and  heat  recovery  efficiency  with 
low  pollution  rate  make  fuel  cells  [1]  potential  useful  energy 
conversion  systems.  Initial  modelling  efforts  have  been 
focused  on  planar  solid  oxide  fuel  cells  (SOFCs).  Experi¬ 
mental  data  are  scarce  and  for  this  reason  substantial  effort  is 
being  devoted  to  developing  numerical  analysis  tools  cap¬ 
able  of  performing  calculations  on  transport  and  electro¬ 
chemical  phenomena  within  the  passages  of  fuel  cells. 

Since  the  first  SOFC  computations,  Vayenas  and  Hegedus 
[2],  the  detail  of  the  mathematical  modeling  has  increased. 
Numerical  simulations  have  been  conducted  at  the  electrode, 
cell,  and  stack  levels.  Modelling  at  the  electrode  level  aims 
at  building  better  electrodes  through  study  of  microscopic 
processes,  while  modelling  at  the  stack  level  aims  at  opti¬ 
mizing  the  design,  by  considering  alternatives  and  determin¬ 
ing  operational  strategies.  Chemical  reactions  (shift  reactions 
and  internal  reforming),  electrical  potential  distribution,  and 
porous-media  flow  are  all  issues  to  be  addressed. 

Fiard  and  Herbin  [3],  Ferguson  [4],  Herbin  et  al.  [5],  and 
Bernier  et  al.  [6]  developed  a  detailed  three-dimensional 
(3D)  SOFC  model  with  governing  equations  for  mass,  heat 
and  electrical  current  for  both  solid  and  gas-channel  flows. 
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Numerical  schemes  for  the  boundary  condition  at  the  inter¬ 
faces  between  the  electrolyte  and  electrodes  are  given.  The  3D 
model  was  applied  to  planar  stack  simulations.  Karoliussen 
and  Nisanciouglu  [7],  Achenbach  [8],  Bessette  and  Wepfer 
[9],  and  Bernier  et  al.  [6]  took  into  account  the  reforming  and 
shift  reactions  in  their  respective  models.  In  addition,  Ahmed 
et  al.  [10],  Sira  and  Ostenstad  [11],  Achenbach  [8],  Bessette 
and  Wepfer  [9],  Costamagna  and  Honegger  [12],  Chan  et  al. 
[13],  Dong  et  al.  [14]  and  Beale  et  al.  [15,16]  applied  their 
models  to  heat  and  mass  transfer  in  SOFCs. 

The  complexity  of  the  SOFC  problem  requires  the  use 
of  large  fast  computers  to  tessellate  the  geometry  into  a 
large  number  of  mesh  points,  and  solve  the  coupled  partial 
differential  equations  describing  the  transport  phenomena. 
The  theoretical  framework  for  stack  modeling  based  on 
simplified  numerical  methods,  so  that  numerical  simulation 
become  tractable  on  personal  computers,  was  introduced  by 
various  authors,  Achenbach  [8],  Bernier  et  al.  [6],  Beale  et  al. 
[15].  Both  single  cells  and  stacks  of  fuel  cells  are  considered 
in  the  present  work.  Fig.  1  is  a  schematic  of  a  stack  con¬ 
sidered  in  this  study. 


2.  Mathematical  formulation 

2.1.  Detailed  numerical  model  (DNM) 

For  single  cells  and  small  stacks  it  is  possible  to  discretize 
the  entire  domain  and  solve  the  governing  equations  directly. 
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This  is  referred  to  below  as  a  detailed  numerical  model 
(DNM).  The  equations  to  be  considered  are  the  usual 
transport  equations,  namely 

=  F §rad 0  +  $  (1) 

where  0  takes  the  value  1  (continuity),  U  (momentum),  yt 
(mass  fraction)  and  h  (enthalpy),  and  r  and  S  are  exchange 
coefficients  and  source  terms,  respectively.  Reynolds  num¬ 
bers  for  both  fuel  and  air  are  small,  and  a  turbulence  model 
was  not  therefore  invoked.  Solid  and  fluid  physical  proper¬ 
ties  were  similar  to  those  given  in  Beale  et  al.  [16]. 

At  the  anode  surface,  electrochemical  oxidation  takes 
place  as 

H2  +  O2-  -f  H20  +  2e"  (2) 

CO  +  O2-  -*•  C02  +  2e“  (3) 

At  the  cathode  surface,  reduction  takes  place 
02  +  4e“  ->  202-  (4) 


The  surface  rates,  J,  for  H2,  H20  and  02  are  required  to 
compute  mass  and  species  source  terms.  These  can  be 
related  to  local  current  density,  i,  by  Faraday’s  law 


J 


M  i 
1000 vF 


(5) 


where  M  is  molecular  weight,  v  the  valence,  and  F  the 
Faraday’s  constant.  The  cell  voltage,  V,  can  be  computed  as 


V  =  E  —  iRj  —  »7a  —  rje  =  E  —  iR / 


(6) 


where  ?/a  and  rjc  are  anodic  and  cathodic  overpotentials.  R; 
(£1  m2)  is  the  local  Ohmic  resistance,  R/  (£2  m2)  can  be 
regarded  as  a  locally  Tumped  internal  resistance’  of  the  cell. 
E  is  the  Nemst  potential; 


£  =  £°4 


yH2y^ 

vh2o 


(7) 


where  y,  are  mole  fractions,  and  Pa  the  air  pressure.  For  the 
results  presented  in  this  paper,  a  fourth-order  least-squares 
polynomial  [14]  fitted  to  experimental  data  in  the  range 
550-1200  °C,  was  used  to  compute  the  Ohmic  resistance 


and  electrolyte  overpotentials  in  Eq.  (6),  though  latterly  a 
Butler-Volmer  equation  is  now  being  used  compute  the 
overpotentials  in  Eq.  (6),  explicitly. 

The  local  heat  source  due  to  the  electrothermal  effect  of 
Ohmic  resistance  and  overpotentials  can  be  expressed  as 


i(E~V) 

Lcell 


(8) 


where  Lcell  is  the  thickness  of  the  cell  sheet. 

A  finite  volume  method  is  used  to  solve  the  partial 
differential  equations  subject  to  appropriate  boundary  con¬ 
ditions.  The  geometry  is,  such  that  a  Cartesian  mesh,  passing 
through  both  solid  and  fluid  zones,  was  conveniently 
employed.  Fig.  2a  shows  a  computational  mesh  used  to 
perform  calculations  using  a  DNM.  The  main  components  of 
the  cell  are  fuel  channels,  electrolyte  and  electrodes,  air 
channels,  and  interconnects.  Fuel  and  air  are  in  cross-flow. 
Two  designs  were  considered:  (a)  with  both  the  fuel  and  air 
channels  in  the  form  of  flat  rectangular  ducts;  (b)  with 
numerous  individual  air  channels  and  a  single  rectangular 
fuel  channel.  Overall  dimensions,  boundary  conditions,  etc. 
are  the  same  in  both  cases.  Current  density,  temperature, 
power  density,  lumped  resistance,  and  fuel  utilization  were 
all  computed. 

The  calculation  proceeds  as  follows:  (1)  initial  values  are 
assumed  for  transport  properties  and  cell  voltage  V.  (2)  The 
main  calculation  procedure  is  commenced  and  heat  and 
mass  source  terms  computed.  The  transport  equations, 
Eq.  (1),  are  solved.  (3)  The  Nernst  potential  and  internal 
resistance  are  then  computed,  and  the  local  current  density 
obtained.  Steps  (2)  and  (3)  are  repeated  until  sufficient 
convergence  is  obtained.  Either  the  cell  voltage,  V,  or  the 
current,  /,  (or  equivalently  mean  current  density,  i )  must  be 
prescribed.  Earlier  work  by  our  group  was  for  the  case  of 
prescribed  cell  voltage,  V.  Our  preference  for  prescribing 
current  density  here  is  driven  by  the  need  to  ensure  charge 
conservation  from  cell-to-cell  within  a  stack.  For  this  case, 
the  cell  voltage  is  adjusted  iteratively  until  the  correct  value 
of  i  is  obtained.  This  is  quite  straightforward. 

The  DNM  was  used  to  calculate  performance  in  both 
single-cells  and  in  the  manifolds  and  passages  of  stacks  of 
cells.  Two  distinct  commercial  computer  codes,  Phoenics 
and  Fluent,  were  employed  for  this  class  of  analysis.  In  both 
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cases  the  details  of  the  electrochemical  reactions  were 
written  in  Fortran  or  C  and  called  during  run-time  by  the 
flow  solver. 

2.2.  Distributed  resistance  analogy  (DRA)  and 
presumed  flow  methods  (PFMs) 

Because  detailed  numerical  simulations  require  very  large 
computational  resources,  an  alternative  methodology  was 
devised  for  stack  modeling.  The  method  is  a  modified 
version  of  the  DRA  of  Patankar  and  Spalding  [17].  The 
flow  of  both  working  fluids  with  the  associated  coupled  heat/ 
mass  transfer  is  computed  using  local  volume-averaging  so 
that, 

+  div(rP3)i  =  ns,  (9) 

=  -n  grad  pi  +  div(V/<  grad  u)i  +  r,S;  (10) 

+  di  v{rp<j>)i  =  di  v(rT  grad  $).  +  (1 1) 

where  i  is  the  air,  fuel,  electrolyte  (including  the  electrodes), 
interconnect  materials,  as  appropriate.  Source  terms  S',  = 
{Fru)i  are  introduced  into  Eq.  (10);  F  is  referred  to  as  a 
‘distributed  resistance’.  For  heat  transfer,  inter-phase  terms 
having  the  form  rtSj  =  —  cf)  are  introduced  to 

account  for  fluid-to-solid  and  solid-to-fluid  heat  transfer. 
Values  of  F  and  a  may  be  obtained  from  analytical,  numerical 
or  experimental  analysis.  Here  analytical  expressions  for  flow 
and  heat  transfer  in  ducts  were  used,  Beale  et  al.  [18]. 

Continuity  and  concentration  sources  and  sinks  per  unit 
volume  account  for  heterogeneous  chemical  reactions.  State 
variables  may  now  be  considered  as  being  interstitial  bulk 
values.  Diffusive  terms  are  eliminated  from  fluid  regions, 
while  convective  terms  are  superfluous  in  solids.  Two  velo¬ 
cities  and  pressures,  corresponding  to  air  and  fuel,  are  solved 
for  in  each  computational  cell,  and  temperatures  in  all  fluid 
and  solid  zones.  This  is  implemented  with  a  multiply  shared 
space  method,  Spalding  and  Zhubrin  [19],  Zhuhrin  [20],  The 
domain  is  decomposed  into  four  separate  spatial  blocks  in 
order  to  solve  for  air,  fuel,  electrolyte,  and  metallic  inter¬ 
connect  materials,  as  shown  in  Fig.  2b.  Inter-phase  source 
terms  are  passed  back  and  forth  spatially  between  the  four 
blocks.  The  computer  code  Phoenics  was  used  to  develop 
this  capability.  At  present  the  DRA  is  only  implemented  for 
constant  current  density. 

A  presumed  flow  method  (PFM),  based  on  rate  equations, 
represents  the  simplest  possible  methodology.  The  pressure- 
corrected  momentum  equations  are  not  solved,  but  the 
continuity  equations,  which  assume  a  simplified  ordinary 
differential  form  for  single-pass  cross-flow,  account  for  the 
production  and  destruction  of  chemical  species.  For  heat/ 
mass  transfer,  inter-fluid  transfer  terms  and  Ohmic  heat  per 
unit  volume  in  the  electrolyte  are  as  prescribed  as  above. 


Beale  et  al.  [16]  showed  that  the  heat  conduction  in 
the  highly-conducting  interconnect  material  cannot  be 
neglected,  and  an  iterative  solver  was  therefore  used  to 
obtain  a  numerical  solution  to  the  governing  system  of 
equations.  Although  referred  to  as  a  ‘presumed  flow’  method 
the  mass  flow  rates  are  in  fact  adjusted  to  account  for  the 
generation/depletion  of  matter  due  to  electrochemical  reac¬ 
tions.  A  rate  equation  is  used,  as  above,  to  compute  wall 
mass  fraction  from  bulk  values  and  hence  wall  mole  frac¬ 
tions  for  use  in  the  Nemst  equation,  Eq.  (7),  and  the  source 
terms  in  the  transport  equations.  A  multiply- shared  method 
is  not  required,  since  an  in-house  PFM  code  was  developed 
in  C++,  where  the  spatial  discretization  is  chosen  to  corre¬ 
spond  to  the  cell  materials.  The  code  is  capable  of  modelling 
single  cells  and  fuel  cell  stacks,  for  variable  current  density, 
however  unlike  the  DRA  method,  flow  in  the  manifolds 
cannot  be  detailed  at  present,  only  that  within  the  passages 
of  the  fuel  cell. 


3.  Results 

All  three  classes  of  code  were  employed  in  the  study  of 
both  single  fuel  cells  and  stacks  of  cells  under  a  variety  of 
operating  conditions.  The  dimensions  of  the  reference  geo¬ 
metry  are  nominally  0.1m  x  0.1  m.  Boundary  conditions 
and  property  values  are  similar  to  those  given  in  reference 
[16].  Detailed  comparisons  of  the  models  were  undertaken 
for  a  single  cell  of  known  geometry  with  constant  mass 
source  (i.e.  current  density)  and  heat  source,  and  also  under 
the  more  realistic  situation  where  average  current  (density) 
is  pre-defined,  and  local  current  density  and  resistance 
computed,  iteratively.  For  a  10-cell  stack,  a  model  for 
variable  local  current  density  is  used  to  perform  calculations 
in  the  absence  of  manifolds,  i.e.  presumed  uniform  flow  at 
fuel  and  air  inlets.  Studies  were  also  conducted  for  full 
manifold- stack  assemblies  corresponding  to  an  actual  design 
based  on  the  assumption  of  constant  local  current  density. 

Fig.  3  shows  level-lines  for  temperature  obtained  from 
DNM  and  PFM  methods  under  the  ‘idealized’  conditions 
of  presumed  constant  mass/heat  source  corresponding  to 


Fig.  3.  Cell  temperature,  I ;  (°C),  constant  i  =  4  000  A/m2,  P  =  1.5  W/m3 
(plan  view). 
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i  =  4  000  A/m2,  q  =  1.5  X  106  W/m3.  Figs.  4-10  show 
results  for  a  single  cell  with  variable  local  current  density 
corresponding  to  a  prescribed  mean  value.  Figs.  4-6  show 
hydrogen  mass  fraction,  cell  temperature,  and  current  den¬ 
sity,  for  i  =  4  000  A/m2;  the  local  cell  resistance  is  com¬ 
puted  according  to  reference  [14].  Data  obtained  from  DNM 
and  PFM  calculations  are  exhibited.  Figs.  7-10  are  similar 
plots  of  cell  temperature,  and  current  density  for  i  =  6  000 
and  2  000  A/m2  (DNM  only).  Fig.  1 1  shows  the  voltage 
versus  current  density  characteristic  curve,  corresponding  to 
H2  and  02  utilisation  factors  of  0.6  and  0.25,  respectively. 

The  results  of  calculations  of  the  performance  of  a  10-cell 
stack  are  shown  in  Figs.  12-15.  Figs.  12  and  13  show 
elevation  views  of  temperature  in  a  10-cell  stack,  assuming 
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constant  i  =  4  000  A/m2,  q  =  1 .5  x  106  W/m3  obtained 
using  the  DNM  and  DRA  methods.  Figs.  14  and  15  are 
views  of  the  temperature  distribution  and  current  density 
obtained  for  variable  current  density  corresponding  to  a 
mean  value  of  i  =  4  000  A/m2.  The  latter  are  for  a  stack 
with  no  manifolds. 


ol - * - ' - ' - 1 
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Fig.  1 1 .  Cell  voltage  as  a  function  of  mean  current  density. 
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Fig.  13.  Temperature,  T,  (°C),  constant  (  =  4000  A/m2,  P=  1.5  W/m3. 
10-cell  stack,  constant  temperature  walls,  Tw  =  750  °C. 


Fig.  15.  Current  density  (A/m2),  i  =  4  000  A/m2,  PFM. 


4.  Discussion 

Inspection  of  Fig.  3  reveals  that  for  the  ‘idealized  case’ 
of  uniform  heating  and  mass  transfer  by  the  electrolyte; 
the  temperature  is  lowest  at  the  location  corresponding  to 
the  (bottom-left)  air-fuel  inlets,  and  highest  at  the  corre¬ 
sponding  outlets  (top-right).  The  bi-linear  temperature 
distribution  associated  with  cross-flow  is  due  to  the  overall 
energy  balance  being  dominated  by  the  convection  and 
heat-source  terms.  Thus,  even  if  the  fluid  flow  and  che¬ 
mical  reaction  rates  are  completely  uniform,  there  will  be 
temperature  gradients  in  the  fuel  cell:  one  factor  tending 
to  smooth  and  reduce  undesirable  gradients  [16]  is  the 
metallic  interconnect.  Thus  for  thermal  management,  a 
thick,  highly-conducting  interconnect  is  desirable.  It  can 
be  seen  from  Fig.  3,  that  there  is  reasonably  good  agree¬ 
ment  between  the  data  for  the  DNM  and  PFM  methodo¬ 
logies.  The  differences  are  due  to  the  approximate  nature 
in  the  choice  of  Nusselt  numbers  required  for  the  PFM 
model. 

Inspection  of  Figs.  4-6,  for  variable  local  current  density 
and  electrical  resistance,  further  reinforce  the  fact  that  a 
reasonable  estimate  of  the  thermo-mechanical  and  electro¬ 
chemical  performance  of  a  fuel  cell  in  cross-flow,  can  be 
obtained  using  a  PFM,  at  a  fraction  of  the  computer  speed 
and  memory  required  to  perform  a  detailed  computational 
fluid  dynamcs  (CFD)  based  calculation,  provided  the  inlet 
conditions  are  uniform.  Contours  of  Nemst  potential  (not 
shown)  follow  those  of  hydrogen  mass  fraction,  Fig.  4. 
Because  current  density,  and  hence  power  density,  are  higher 
at  the  top-left  comer  of  Fig.  6,  the  global  temperature 
maximum  shifts  to  this  location.  The  cell  resistance  is  also 
a  minimum  at  this  point.  Non-uniform  reaction  kinetics 
will  re-distribute  the  temperature  field,  but  may  not  neces¬ 
sarily  result  in  higher  temperature  gradients  (NB:  the  results 
of  Figs.  3  and  5  are  not  quantitively  comparable,  due  to 
differences  in  cell  resistance).  Reference  [18]  discusses  how 
uniform  flow  conditions  may  be  achieved.  At  i  =  4  000  A/m2, 
the  cell  voltage  computed  from  the  PFM  is  within  0.5%  of 
that  obtained  with  the  DNM,  with  a  lower  cell  voltage 
predicted  by  the  DNM,  and  corresponding  heat  source 
in  Eq.  (8),  thus  affecting  the  temperature  distribution  in 
Fig.  3,  though  agreement  between  DNM  and  PFM  is  still 
within  5  °C. 

Detailed  calculations  were  also  performed  for  average 
current  densities  of  2  000  and  6  000  A/m2.  Agreement  in 
terms  of  local  current  density,  between  the  PFM  and  DNM 
distribution  (not  shown)  is  generally  good.  Comparison  of 
the  results  for  i  =  2  000,  4  000,  6  000  reveals  the  tempera¬ 
ture  distribution  to  be  a  strong  function  of  the  load,  while  the 
current  density  is  relatively  independent,  indicating  that  low 
fuel  and  air  utilizations  are  achieved.  At  i  =  2  000  A/m2, 
the  temperature  profile  resembles  that  of  a  simple  cross- 
flow  heat  exchanger,  due  to  the  difference  in  temperature 
between  the  two  working  fluids  being  large  compared  to 
the  overall  temperature  rise  due  to  the  electrical  heating. 
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At  i  =  6  000  A/m2,  the  converse  is  true.  The  current  density 
contours,  as  shown  in  Figs.  6,8,10  are  relatively  similar  over 
this  range  of  loads.  For  the  DNM  calculations,  very  good 
agreement  was  obtained  when  using  two  different  CFD 
codes;  Phoenics  and  Fluent. 

Inspection  of  the  results  of  the  voltage  current-density 
distribution  curve,  Fig.  11,  also  reveal  that  very  good 
agreement  is  obtained  between  the  DNM  and  PFM  methods, 
at  least  for  values  of  current  density  up  to  i  =  10  000  A/m2. 
At  higher  values  of  i  some  significant  error  becomes  appar¬ 
ent,  due  to  the  assumption  of  constant  Nusselt/Sherwood 
number  being  invalid  at  these  high  mass  transfer  rates. 

Figs.  12  and  13  show  temperature  distributions  for  a  10- 
cell  stack  model  at  constant  i  =  4  000  A/nr  under  adiabatic 
and  isothermal  boundary  conditions,  respectively.  The 
results  shown  in  Fig.  12  are  for  presumed  adiabatic  walls, 
corresponding  to  a  well-insulated  stack.  Excellent  agree¬ 
ment  is  observed  between  the  results  of  the  DNM  and  DRA 
calculations.  The  ‘zig-zag’  pattern  exhibited  by  the  DNM  is 
absent  for  the  DRA  because  a  multi-space  approach  is 
employed,  and  only  the  air-space  temperature  is  displayed, 
whereas  for  the  DNM,  the  temperatures  in  the  air,  fuel, 
electrolyte  and  interconnect  zones  vary  on  a  cell-by-cell 
basis.  The  arising  temperature  field  is  three-dimensional 
(3D)  [16]  even  though  the  fluid  flow  and  current  density 
distributions  are  uniform.  This  is  due  to  the  fact  that  the  fuel 
and  air  streams  are  at  different  temperatures,  and  that  the 
order  of  the  materials  is  repeated  in  the  vertical  direction. 
Thus  the  common  assumption  that  a  single-cell  model 
with  periodic  boundary  conditions  may  be  used  to  predict 
stack  performance  is  erroneous.  In  the  current  DRA  imple¬ 
mentation,  the  vertical-k-direction  influence  was  correctly 
accounted-for  by  meshing  the  geometry  so  that  computa¬ 
tional  cells  corresponded  to  fuel  cells,  and  prescribing  the 
electrode-fuel  (ef)  pair  of  inter-phase  source  terms  in 
Eq.  (11)  so  that  these  were  across  n  —  1  neighbour  values 
aef(0e,ifc+i  —  0r,/t)  and  aef(</>f  —  </>e  fc),  not  n  in-cell  values 
D'etre  k  —  0f  k)’  k  =  1,2, 3 as  were  all  other  terms 
(electrode-air,  air-interconnect,  fuel-interconnect,  etc.). 
With  this  important  modification,  it  can  be  seen  that  the 
DRA  approach  generates  the  3D  results  obtained  with  the 
detailed  CFD  simulation  at  a  fraction  of  the  computational 
cost.  Any  volume-averaging  or  ‘porous-media’  technique 
which  does  not  account  for  the  effects  due  to  the  ordering  of 
the  streams,  will  incorrectly  generate  constant  temperature 
profiles  in  the  vertical  direction. 

Fig.  13  shows  the  results  for  the  corresponding  constant 
wall  temperature  case,  corresponding  to  the  limiting  case  of 
an  uninsulated  stack  immersed  inside  an  oven  at  750  °C. 
The  reader  will  note  that  although  the  isotherms  at  the  wall 
in  Fig.  13  are  parallel  rather  than  normal  to  the  wall,  similar 
temperature  distributions  are  observed  in  the  interior  of 
the  fuel  cell  as  for  the  adiabatic  case,  Fig.  12.  This  suggests 
that  temperature  control  may  be  a  matter  for  concern  for 
the  fuel  cell  designer  employing  prototypes  in  experimen¬ 
tal  test  rigs,  where  the  interior  temperatures  may  exceed 


surface  values.  Figs.  14  and  15  show  the  contours  of 
temperature  and  current  density  in  a  10-cell  stack  obtained 
at  variable  current  density,  corresponding  to  i  =  4  000  A/m2. 
The  secondary  temperature  effects  are  less  pronounced 
than  in  Fig.  13,  because,  the  power  consumption  is  lower; 
however,  they  are  readily  apparent.  The  reader  will  also 
note  that  predicted  local  current  density  displays  3D 
effects,  due  to  the  sensitivity  of  cell  resistance  to  local 
temperature. 


5.  Conclusions 

Calculations  were  performed  on  single  and  10-cell  stacks 
of  SOFC  fuel  cells  using  three  distinct  approaches  referred 
to  as  the  PFM,  DRA  and  DNM.  Both  constant  heat  and 
current  density,  and  variable  local  current  density,  corre¬ 
sponding  to  known  average  values  of  i,  were  considered 
under  adiabatic  and  isothermal  wall  conditions.  For  single 
cells,  all  of  these  methods  can  be  reliably  used  to  perform 
calculations  in  planar  SOFCs  with  the  DNM  being  the 
most  accurate  since  it  does  not  require  prescription  of 
overall  heat  and  mass  transfer  coefficients.  For  large 
stacks  of  fuel  cells,  however  the  DNM  requires  very  large 
computational  meshes,  resulting  in  large  data  sets  and 
compute  times.  Under  these  circumstances,  if  it  is  known 
that  the  flow  is  fully-developed  and  uniform,  a  PFM  should 
be  considered  as  an  alternative  or  complementary  simula¬ 
tion  tool  at  a  fraction  of  the  compute  cost.  The  DRA 
represents  a  compromise  between  cost  and  performance, 
and  can  be  used  in  situations  where  the  flow  and  pressure 
distributions  in  the  stack  are  not  necessarily  uniform.  It  is 
important  that  when  ‘local  volume-averaging’  techniques 
are  employed,  secondary  transport  phenomena  are  not 
obscured. 

The  key  difference  between  the  PFM/DRA  and  the  DNM 
program  is  in  the  specification  of  overall  heat/mass  coeffi¬ 
cients.  Any  errors  in  these  values  will  impact  the  Nernst 
equation  and  the  local  cell  resistance.  The  results  presented 
in  this  paper  are  relatively  preliminary,  and  we  anticipate 
further  reconciliation  between  the  different  modelling 
techniques  as  experience  is  gained.  It  is  to  be  noted  that 
many  alternatives  are  possible;  the  DRA  could  be  replaced 
by  a  code  whereby  detailed  flow  models  based  on  CFD 
are  obtained  in  necessary  zones,  such  as  manifolds,  but 
that  a  PFM  model  be  adopted  in  the  stack  where  addi¬ 
tional  storage  is  provided-for  auxiliary  variables,  such 
as  current  density  and  potential,  which  are  not  required 
elsewhere. 


6.  Future  work 

A  semi-empirical  resistance  model,  which  combines  the 
Ohmic  and  overpotential  terms  was  used  in  this  paper,  owing 
to  the  sparse  experimental  data  available  for  the  particular 
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design  under  consideration.  This  has  now  been  supplanted 
with  a  Butler- Volmer  equation  for  anodic  and  cathodic  over¬ 
potential  in  most  of  the  codes,  and  the  latter  will  be  adopted 
from  now  on.  To  further  improve  the  correlation  between 
the  PFM,  DRA,  and  DNMs,  additional  effort  is  required 
when  prescribing  mass  transfer  coefficients,  and  use  of 
Stefan-Maxwell  equations  for  multi-species  mass  diffusion 
coefficients,  in  place  of  the  simple  Fickean  analysis  deve¬ 
loped  here,  is  a  worthwhile  undertaking. 

Many  workers  are  now  modelling  the  electrical  potential 
and  current  collection  associated  with  porous-media  diffu¬ 
sion  and  chemical  reactions  in  a  fully-coupled  manner.  The 
simplified  model  used  here,  where  the  electrolyte  includes 
the  electrodes  is  considered  appropriate  for  thin  planar 
geometries,  but  may  not  be  valid  for  other,  more  complex 
fuel-cell  designs.  Under  these  circumstances  detailed 
numerical  analysis  may  subsequently  be  used  to  prescribe 
electrical  and  thermal  shape  factors,  needed  to  compute 
heat/mass  transfer  coefficients  for  the  rate  equations.  More¬ 
over  the  details  of  the  chemical  reactions  have  been  kept 
simple;  for  a  first  analysis  important  phenomena  such  as 
internal  reforming  have  been  excluded  in  the  interest  of 
simplicity.  For  the  same  reasons  radiation  heat  transfer  was 
also  neglected.  Future  work  will  include  these. 

The  success  of  the  PFM  and  the  DRA  are  due,  in  part  to 
the  ease  of  constructing  models  for  single-pass  cross-flow. 
For  multi-pass  flow  and/or  where  the  flow  channels  are  of 
arbitrary  location  and  direction  (i.e.  not  necessarily  oriented 
along  well-defined  geometric  lines),  the  scientific  principles 
remain  the  same,  but  the  problem  of  engineering  codes 
under  these  circumstance  is  a  challenge.  In  modifying  the 
DRA  to  correctly  predict  3D  phenomena,  it  was  necessary 
for  the  computational  mesh  to  conform  to  the  boundaries  of 
the  fuel  cells  in  the  stack,  so  that  neighbour  values  be 
correctly  prescribed.  Some  further  research  work  is  required 
to  allow  for  the  DRA  to  be  employed  with  aribitrary  non¬ 
body-fitting  meshes. 

The  need  for  reliable  experimental  data  is  readily  appar¬ 
ent,  and  while  gathering  such  data  can  be  substantially  more 
time-consuming  than  constructing  models  and  performing 
numerical  calculations,  the  availability  of  such  data  are 
critical  for  the  evaluation  and  improvement  of  existing 
numerical  models.  Although  most  of  the  work  to-date  has 
been  concerned  with  SOFCs,  much  of  the  modelling  exper¬ 
tise  can  readily  be  modified  for  application  to  other  fuel 
cells,  such  as  proton  exchange  membrane  fuel  cells.  This 
avenue  is  actively  being  pursued  by  a  majority  of  the 
authors. 
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